## POQ Replication for Trends Analyses
## by: AR Flores - for Lewis et al.

## CONTACT

df1 <- haven::read_dta("C:/Users/aflor/Downloads/Personal Contact Time Series.dta")

df1 <- df1 %>%
  filter(questiontype!="contact")


df1 <- df1 %>%
  mutate(s_star = round((contact/100)*sample, 0),
         date_df = date,
         date = as.numeric(date))

df1 <- df1 %>%
  mutate(int_method = case_when(survey=="PRRI" | survey=="PRRI/Princeton" | survey=="CNN/ORC" ~ 0,
                                survey=="KU" | survey=="KU/SSI" | survey=="Pew" | survey=="Williams/IPSOS" ~ 1),
         wording = case_when(survey =="CNN/ORC" | survey=="KU" ~ "Do you happen to have a family member or close friend who is transgender?",
                             survey=="KU/SSI" | survey=="PRRI/Princeton" ~ "Do you have a close friend of family member who is transgender?",
                             survey=="Pew" ~ "Thinking about the people or people you know who are transgender, are they a close friend or family member?",
                             survey=="PRRI" ~ "Now thinking about the people that you know. Please tell me whether you have a close friend or family members who is trangender.",
                             survey=="Williams/IPSOS" ~ "How familiar are you with transgender people?"                             ))

md1 <- bf(s_star | trials(sample) ~ s(date, k=6) + int_method + (1|wording))

df1 <- as.data.frame(df1)

ft1 <- brm(md1, data = df1, family = binomial(link="logit"), warmup=1000, iter=2000)

ft1_1 <- brm(md1, data = df1, family = binomial(link="logit"), warmup=1000, iter=2000, chains=4, cores = 4)

ft1_2 <- brm(md1, data = df1, family = binomial(link="logit"), warmup=1000, iter=2000, chains=4, cores = 4,backend = "cmdstanr")


summary(ft1)

sm2 <- conditional_effects(ft1)

sm2$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")

p <- ggplot()

p + geom_ribbon(data=sm2$date, aes(x = newdate, ymin=lower__, ymax = upper__), alpha=0.2) +
  geom_point(data=df1, aes(x=date_df, y = contact/100), size=2, shape=1) +
  geom_line(data=sm2$date, aes(x = newdate, y = estimate__), size = 1.2) +
  scale_y_continuous(labels=scales::percent, limits = c(0,0.25), breaks=seq(0,.25,.05)) +
  scale_x_date() +
  ggthemes::theme_clean() +
  labs(x = NULL,
       y = "Percent")


df1$lb <- (df1$contact/100) - qnorm(.975)*sqrt(((df1$contact/100)*(1-df1$contact/100))/df1$sample)
df1$ub <- (df1$contact/100) + qnorm(.975)*sqrt(((df1$contact/100)*(1-df1$contact/100))/df1$sample)

p <- ggplot()

p + geom_ribbon(data=sm2$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray") +
  geom_line(data=sm2$date, aes(x = newdate, y = estimate__), size = 1.2, color="red") +
    geom_pointrange(data = df1, aes(x = date_df, y = contact/100, ymin = lb, ymax=ub), color="black") +
  scale_y_continuous(labels=scales::percent, limits = c(0,0.35), breaks=seq(0,.25,.05)) +
  scale_x_date() +
  ggthemes::theme_clean() +
  labs(x = NULL,
       y = "Percent")

ggsave("contact_controls.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", dpi=300, width=4, units="in")
## Military

df2 <- haven::read_dta("C:/Users/aflor/Downloads/Military Time Series 2015-21.dta")

#this one does all levels / 3 models.
df2 <- df2 %>%
  mutate(s_allow = round((allow/100)*sample,0),
         s_not = round((notallow/100)*sample,0),
         s_dk = round((dk/100)*sample,0),
         date_df = date,
         date = as.numeric(date))
df2 <- df2 %>%
  mutate(int_method = case_when(var2=="Gallup" | var2=="UI-Springfield" | var2=="Quinnipiac" | var2=="UD" | var2=="PRRI" ~ 0,
                                var2=="KU" | var2=="KU/SSI" | var2=="GfK" | var2=="Pew" ~ 1),
         wording = case_when(var2=="KU" | var2=="GfK" | var2=="KU/SSI" ~ "Transgender people should be allowed to serve openly in the military.",
                             var2=="UI-Springfield" | var2=="UD" | var2=="PRRI" ~ "Do you support or oppose allowing transgender persons to serve openly in the US military?",
                             var2=="Quinnipiac" ~ "Do you think transgender people should be allowed to serve in the military or not?",
                             var2=="Pew" ~ "Banning transgender people from serving in the military.",
                             var2=="Gallup" ~ "Do you favor or oppose transgender men and women to serve in the military?"
                             ))

md2 <- bf(s_allow | trials(sample) ~ s(date, k=6)+ int_method + (1|wording))
md3 <- bf(s_not | trials(sample) ~ s(date, k=6)+ int_method + (1|wording))
md4 <- bf(s_dk | trials(sample) ~ s(date, k=6)+ int_method + (1|wording))

df2 <- as.data.frame(df2)

ft2 <- brm(md2, data = df2, family = binomial(link="logit"), warmup=1000, iter=2000)
ft3 <- brm(md3, data = df2, family = binomial(link="logit"), warmup=1000, iter=2000)
ft4 <- brm(md4, data = df2, family = binomial(link="logit"), warmup=1000, iter=2000)

summary(ft2)
summary(ft3)
summary(ft4)

sm2 <- conditional_effects(ft2)
sm3 <- conditional_effects(ft3)
sm4 <- conditional_effects(ft4)

sm2$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")
sm3$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")
sm4$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")

df2$lb_al <- (df2$allow/100) - qnorm(.975)*sqrt(((df2$allow/100)*(1-df2$allow/100))/df2$sample)
df2$ub_al <- (df2$allow/100) + qnorm(.975)*sqrt(((df2$allow/100)*(1-df2$allow/100))/df2$sample)

df2$lb_n <- (df2$notallow/100) - qnorm(.975)*sqrt(((df2$notallow/100)*(1-df2$notallow/100))/df2$sample)
df2$ub_n <- (df2$notallow/100) + qnorm(.975)*sqrt(((df2$notallow/100)*(1-df2$notallow/100))/df2$sample)

df2$lb_d <- (df2$dk/100) - qnorm(.975)*sqrt(((df2$dk/100)*(1-df2$dk/100))/df2$sample)
df2$ub_d <- (df2$dk/100) + qnorm(.975)*sqrt(((df2$dk/100)*(1-df2$dk/100))/df2$sample)

p <- ggplot()

p + geom_ribbon(data=sm2$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray") +
  geom_ribbon(data=sm3$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray", alpha=0.6) +
  geom_ribbon(data=sm4$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray", alpha=0.6) +
  geom_line(data=sm2$date, aes(x = newdate, y = estimate__, linetype = "Allow"), size = 1, color="black") +
  geom_line(data=sm3$date, aes(x = newdate, y = estimate__, linetype = "Not Allow"), size = 1, color="red") +
  geom_line(data=sm4$date, aes(x = newdate, y = estimate__, linetype = "DK/Refused/Skipped"), size = 1, color="blue") +
  geom_pointrange(data = df2, aes(x = date_df, y = allow/100, ymin = lb_al, ymax=ub_al, shape="Allow"), color="black", size = 1, fatten = 1.5) +
  geom_pointrange(data = df2, aes(x = date_df, y = notallow/100, ymin = lb_n, ymax=ub_n, shape = "Not Allow"), color="red", size=1, fatten = 1.5) +
  geom_pointrange(data = df2, aes(x = date_df, y = dk/100, ymin = lb_d, ymax=ub_d, shape = "DK/Refused/Skipped"), color="blue",  size=1, fatten=1) +
  scale_y_continuous(labels=scales::percent, limits = c(0,.8), breaks=seq(0,.8,.1)) +
  scale_x_date() +
  scale_linetype_discrete("") +
  scale_shape_discrete("") +
  ggthemes::theme_clean() +
  theme(legend.background = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 8)) +
  guides(shape = guide_legend(override.aes = list(size = 0.4, fatten = 1, color=c("black","blue","red")))) +
  labs(x = NULL,
       y = "Percent")

ggsave("military.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", dpi=300, width=4, units="in")


## Thermometers

df3 <- haven::read_dta("C:/Users/aflor/Downloads/Thermometer Time Series.dta")

# weight the regression? 

df3 <- df3 %>%
  mutate(wt = sample / mean(sample),
    date_df = date,
         date = as.numeric(date))

df3 <- df3 %>%
  mutate(int_method = case_when(org == "HRC/Greenberg Quinlan" | org=="ANES" | org=="ANES Post Election" ~ 0,
                                org == "KU" | org=="ANES Pilot/YouGov" | org=="GfK" | org=="KU/SSI"  ~ 1),
         wording = case_when(org=="HRC/Greenberg Quinlan" & date==11887 ~ "a",
                             org=="HRC/Greenberg Quinlan" & date!=11887 ~ "b",
                             org=="GfK" | org=="KU" | org=="KU/SSI" ~ "c",
                             org=="ANES Pilot/YouGov" | org=="ANES Post Election" | org=="ANES" ~ "d"))

md5 <- bf(temp | weights(wt) ~ s(date, k=6) + int_method + (1|wording))

df3 <- as.data.frame(df3)

ft5 <- brm(md5, data = df3, warmup=1000, iter=2000)

summary(ft5)

sm5 <- conditional_effects(ft5)

sm5$date$newdate = as.Date(sm5$date$effect1__, origin = "1970-01-01")

df3$se <- 38.5/sqrt(df3$sample)

df3$lb <- df3$temp - qt(.975, df = df3$sample-1)*df3$se
df3$ub <- df3$temp + qt(.975, df = df3$sample-1)*df3$se

p <- ggplot()

p + geom_ribbon(data=sm5$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray") +
  geom_line(data=sm5$date, aes(x = newdate, y = estimate__), size = 1.2, color="red") +
  geom_pointrange(data = df3, aes(x = date_df, y = temp, ymin = lb, ymax=ub), color="black") +
  scale_y_continuous(limits = c(20,80), breaks=seq(20,80,10)) +
  scale_x_date() +
  ggthemes::theme_clean() +
  labs(x = NULL,
       y = "Mean")

ggsave("therm_controls.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", dpi=300, width=4, units="in")

## Bathrooms

df4 <- haven::read_dta("C:/Users/aflor/Downloads/Bathroom Time Series.dta")

df4 <- df4 %>%
  mutate(s_s = round((supportGI/100)*var8,0),
         s_o = round((opposeGI/100)*var8,0),
         s_dk = round((neither/100)*var8,0),
         date_df = date,
         date = as.numeric(date))

df4 <- df4 %>%
  mutate(int_method = case_when(survey=="CNN/ORC"| survey=="Gallup" | survey=="PT/Marist" | survey=="CBS/NYT" | survey=="NBC/WSJ" | survey=="AP/NORC" | survey=="UD" | survey=="PRRI" | survey=="ANES"| survey=="ANES Pre Election" | survey=="Pew" ~ 0,
                                survey=="KU" | survey=="Williams/IPSOS" | survey=="KU/SSI" | survey=="GfK" ~ 1),
         wording = case_when(survey=="KU" | survey=="GfK" | survey=="KU/SSI" ~ "Transgender people should only be allowed to use public restrooms that are consistent with the sex listed on their driver's license/state ID card",
                             survey=="CNN/ORC" ~ "Overall, would you say you favor or oppose laws that require transgender individuals to use facilities that correspond to their gender at birth rather than their gender identity?",
                             survey=="Gallup" ~ "In terms of policies governing public restrooms, do you think these policies should -- [ROTATED: require transgender individuals to use the restroom that corresponds with their birth gender (or should these policies) allow transgender individuals to use the restroom that corresponds with their gender identity]?",
                             survey=="PT/Marist" ~ "Do you think transgender or gender fluid people should be allowed to choose the public restroom with which they identify or should be required to follow the sex on their birth certificate in using a public restroom?",
                             survey=="CBS/NYT" | survey=="AP/NORC" | survey=="ANES Pre Election" ~ "Do you think people who are transgender--that is, someone who identifies themselves as the sex or gender different from the one they were born as--should be allowed to use the public bathrooms of the gender they identify with or should they have to use the public bathrooms of the gender they were born as?",
                             survey=="NBC/WSJ" ~ "Do you believe transgender people should be allowed to use the public restroom of the gender with which they identify or should they be legally prevented from doing so or do you not have an opinion one way or the other?",
                             survey=="UD" ~ "Laws that require people to use restrooms that match the sex listed on their birth certificate.",
                             survey=="Williams/IPSOS" ~ "Transgender people should  be allowed to use the restroom associated with their gender identity.",
                             survey=="PRRI" ~ "Laws that require transgender individuals to use bathrooms that correspond to their sex at birth rather than their current gender identity.",
                             survey=="Pew" ~ "Transgender people should be.[Allowed to use the public restrooms of the gender with which they currently identify] or [Required to use the public restrooms of the gender they were born into]",
                             survey=="ANES" ~ "Should transgender people... have to use the bathrooms of the gender they were born as, or should they be allowed to use the bathrooms of their identified gender?"))

md6 <- bf(s_s | trials(var8) ~ s(date, k=6) + int_method + (1|wording))
md7 <- bf(s_o | trials(var8) ~ s(date, k=6)+ int_method + (1|wording))
md8 <- bf(s_dk | trials(var8) ~ s(date, k=6)+ int_method + (1|wording))

df4 <- as.data.frame(df4)

ft2 <- brm(md6, data = df4, family = binomial(link="logit"), warmup=1000, iter=2000)
ft3 <- brm(md7, data = df4, family = binomial(link="logit"), warmup=1000, iter=2000)
ft4 <- brm(md8, data = df4, family = binomial(link="logit"), warmup=1000, iter=2000)

summary(ft2)
summary(ft3)
summary(ft4)

sm2 <- conditional_effects(ft2)
sm3 <- conditional_effects(ft3)
sm4 <- conditional_effects(ft4)

sm2$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")
sm3$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")
sm4$date$newdate = as.Date(sm2$date$effect1__, origin = "1970-01-01")

df4$lb_al <- (df4$supportGI/100) - qnorm(.975)*sqrt(((df4$supportGI/100)*(1-df4$supportGI/100))/df4$var8)
df4$ub_al <- (df4$supportGI/100) + qnorm(.975)*sqrt(((df4$supportGI/100)*(1-df4$supportGI/100))/df4$var8)

df4$lb_n <- (df4$opposeGI/100) - qnorm(.975)*sqrt(((df4$opposeGI/100)*(1-df4$opposeGI/100))/df4$var8)
df4$ub_n <- (df4$opposeGI/100) + qnorm(.975)*sqrt(((df4$opposeGI/100)*(1-df4$opposeGI/100))/df4$var8)

df4$lb_d <- (df4$neither/100) - qnorm(.975)*sqrt(((df4$neither/100)*(1-df4$neither/100))/df4$var8)
df4$ub_d <- (df4$neither/100) + qnorm(.975)*sqrt(((df4$neither/100)*(1-df4$neither/100))/df4$var8)

p <- ggplot()

p + geom_ribbon(data=sm2$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray", alpha=0.5) +
  geom_ribbon(data=sm3$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray", alpha=0.5) +
  geom_ribbon(data=sm4$date, aes(x = newdate, ymin=lower__, ymax = upper__), fill="light gray") +
  geom_line(data=sm2$date, aes(x = newdate, y = estimate__, linetype = "Current Gender Identity"), size = 1, color="black") +
  geom_line(data=sm3$date, aes(x = newdate, y = estimate__, linetype = "Birth Gender"), size = 1, color="red") +
  geom_line(data=sm4$date, aes(x = newdate, y = estimate__, linetype = "DK/Refused/Skipped"), size = 1, color="blue") +
  geom_pointrange(data = df4, aes(x = date_df, y = supportGI/100, ymin = lb_al, ymax=ub_al, shape="Current Gender Identity"), color="black", size = 1, fatten = 1.5) +
  geom_pointrange(data = df4, aes(x = date_df, y = opposeGI/100, ymin = lb_n, ymax=ub_n, shape = "Birth Gender"), color="red", size=1, fatten = 1.5) +
  geom_pointrange(data = df4, aes(x = date_df, y = neither/100, ymin = lb_d, ymax=ub_d, shape = "DK/Refused/Skipped"), color="blue",  size=1, fatten=1) +
  scale_y_continuous(labels=scales::percent, limits = c(0,.8), breaks=seq(0,.8,.1)) +
  scale_x_date() +
  scale_linetype_discrete("") +
  scale_shape_discrete("") +
  ggthemes::theme_clean() +
  theme(legend.background = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 8)) +
  guides(shape = guide_legend(override.aes = list(size = 0.4, fatten = 1, color=c("red", "black", "blue")), nrow=2),
         linetype = guide_legend(nrow=2)) +
  labs(x = NULL,
       y = "Percent")

ggsave("bathroom_controls.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", dpi=300, width=4, units="in")
